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Abstract 

Global asymptotic stability of rational difference equations is an area of research that has been well 
studied. In contrast to the many current methods for proving global asymptotic stability, we propose 
an algorithmic approach. The algorithm we summarize here employs the idea of contractions. Given a 
particular rational difference equation, defined by a function Q : R'^^^ R^+^, we attempt to find a 
K value for which shrinks distances to the difference equation's equilibrium point. We state some 
general results that our algorithm has been able to prove, and also mention the implementation of our 
algorithm using Maple. 

1 Introduction 

In this paper we will introduce an algorithmic approach to proving global asymptotic stability (GAS) 
of equilibrium points of rational difference equations. The field of rational difference equations has 
applications to many other fields including biology, economics, and dynamical systems. In application 
areas, one often studies time-evolving sequences produced by recurrences with the goal of discovering end 
behavior of the sequence, given some initial conditions. There are many types of end behavior that one 
may be interested in. We will be concerned only with global asymptotic stability. Essentially, given a 
fixed rational difference equation, when the sequence that it produces converges for any reasonable initial 
conditions, we say that it is GAS. In Section pTT] we will state the precise definition for GAS, as well 
as introduce all of the other definitions necessary to study stability of difference equations. We will also 
state a theorem, originally proved in 8 , which will be the basis of our algorithm. The main algorithm 
is presented in Sections [2] and [3] In the algorithm we first reduce the problem of GAS to the problem 
of proving that a particular polynomial is positive. Then, we prove that a multivariate polynomial is 
positive (when all of its variables are taken to be positive) using our new algorithm. Next, Section [4] 
contains a proof-of-concept that our algorithm is indeed applicable to prove GAS. In Section |5] we state 
a few of the results that our algorithm can prove. In addition, in Section |6] we mention the most useful 
commands in the Maple package that accompanies this paper. 

1.1 Definitions 

Following the various works of Ladas, et. al. [2|[71|9], we begin by stating a few standard definitions 
needed to study rational difference equations and stability. 

Definition 1.1. A rational difference equation (of order /c + 1) is an equation of the form 

Xn-\-l — R{Xri ^Xfi — ■fXn — k') (l) 

where the function R{uo, ui, . . . ,Uk) is a rational function which maps /^+^ to /, for some interval / C R. 
Typically, we will take / to be [0, oo) or (0, oo). 

Given a function R we say that a solution of ([T]) is a sequence {xn}^=_k which satisfies ([T]). One can 
also think of a solution, {a:n}^_/., as being associated to the specific initial conditions . . . ,xo} 

created by repeatedly applying R. If a solution is constant, Xn = x, for all n > —k then we say that the 
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solution is an equilibrium solution, and x is called an equilibrium point, or simply an equilibrium of F. In 
practice, we find the equilibria by solving the equation x = R{x, . . . , x), and taking the solutions which 
lie in the interval /. 

The main topic to be investigated in this paper is end behavior, specifically stability, of a solution of 
a given difference equation. There are various notions of stability that will now be defined. 

Definition 1.2. An equilibrium point, x, of ([T]) is said to be 

1. locally stable if for every e > there exists ^ > such that if {xn}^^_f. is a solution to ([T]) with 
the property that 

\x-k - x| + - x| H \- \xo - x\ < 6 

then \xn — x\ < e for all n > 0. 

2. locally asymptotically stable (LAS) if x is locally stable, and if there exists a 7 > such that if 
{xn}^^_j^ is a solution to ([T]) with the property that 

\x-k - x| + - x| H \-\xo - x\ < 7 

then 

lim Xn = X 

3. a global attractor if for every solution, {xn}^_fc, of ([T]) we have 

lim Xn = X 

n— )-oo 



4. globally asymptotically stable (GAS) if x is a global attractor, and x is locally stable. 

5. unstable if x is not locally stable. 

Our goal in this paper is to present an algorithm to prove GAS. Since GAS implies LAS, the first 
step must be to prove LAS (since, if a difference equation is not LAS it can't be GAS). The linearized 
stability theorem, which provides easily verifiable criteria for local asymptotic stability, can be found in 
many books and papers [3l [4j |9l [TT]. Because it is not central to our algorithm, we will omit the 
theorem and notation needed to state it. 

In contrast to local asymptotic stability which is relatively easy to verify using the linearized stability 
theorem, global asymptotic stability has no similarly general necessary and sufficient conditions. There 
are a handful of theorems, providing sufficient conditions, that have been used to verify the global 
asymptotic stability of many specific difference equations. However, given a difference equation defined 
by the function R, it is not always obvious which theorem to apply. For a discussion of many of these 
theorems see [2 . 

Our algorithm will only rely on the following theorem which is first presented in a paper by Kruse and 
Nesemann 0. It will be stated it in a slightly different manner than it appears in their paper, using the 
notation we have established in this paper. First, it will be necessary to consider the difference equation 
associated to a function R in vector form. Let Q : /^+^ jk+i defined from R as 
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(2) 



Note that this transformation from R to Q essentially creates an order 1 mapping out of an order k -\- 1 
mapping. In addition, Q is now a map that can be composed with itself, so 

where Ao = (xo, . . . , x_fc) is the vector of initial conditions. Now we can state the theorem. 

Theorem 1.1 (Kruse, Nesemann 1999). Let ||-|| denote the Euclidean norm (i.e., ||(a,6)|| - 
Let S denote either [0, cxo) or (0, cxd) (the function Q will necessitate which). Let Q : 'B^'^'^ 



Va^ + b^). 

> 8^^+^ be a 
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continuous mapping of the form Q with a unique fixed point X G S^^"*^. Suppose for the discrete dynamic 
system 

Xn+i = Q{Xn), n = 0,1,2,... (3) 
there exists an integer K > 1 such that the K^^ iterate of Q satisfies 

- A'll < \\X - X\\ for all X G 7^ X. (4) 

Then X is GAS with respect to the norm ||-||. 

First, notice that this integer K tells us which power of Q is a contraction with respect to X, i.e., 
shrinks distances to X. This gives an intuitive reason for ||Q^(A') — X\\ < \\X — X\\ to imply global 
asymptotic stability. The proof of this theorem can be found in [8]. 

Notice that the various definitions of stability, as they were stated in Definition |1.2| do not quite 
apply here because our unique fixed point (or equilibrium) is a vector rather than a scalar. However, 
Definition 1.2 can be easily translated to the vector case. The recurrence is ([3|, the equilibrium is a 



vector solution to the equation Q{X) = X, and the order of the recurrence, /c + 1, is 1 (so /c = 0). Other 



than these minor changes, a word for word translation of Definition 1.2 is what we mean by X being 
GAS in Theorem O 



Next we will see how we utilized this theorem to create a global asymptotic stability proof algorithm. 



2 From Global Asymptotic Stability to Polynomial Posi- 
tivity 

In this section we will see how to reduce the question of global asymptotic stability of a rational difference 
equation to a question about an associated polynomial being positive. Throughout this section assume 
that we have fixed a rational difference equation. 



R{Xn-, • • • 5 Xn — k)i (5) 



of order /c + 1, with a unique equilibrium x. Also assume that is a rational function with positive 
coefficients, so R : [O^oo)^^^ [0,oo), and x is non-negative (if there is no constant term in the 
denominator of R we cannot allow to be in the domain, so R : (0,00)^^^^ (0,oo), and x must be 
strictly positive). In order to apply Theorem |1.1| we must think of ([s} and its equilibrium in their vector 
forms. For example, if Xn+i — /^^^^ then 
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The goal will be to find a positive integer, K, which satisfies Q. Motivated by this goal, we will 
construct the following polynomial, given specific Q, X, and K (assume we have conjectured some value 
for K): 

Pq,x,k{^) = numerator (^\\X - xf - \\q^ (X) - . (6) 

Consider the implication of Pq > ior X > (or > 0, both componentwise), and X X. 

< numerator (^\\X - X\\^ - ||q^(A') - 

^ < llA' - x\\^ - \\q^{x) - 

|b^(A') - A'll < (7) 
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Of course, the first implication, undoing the numerator from line 1 to line 2, in general will not preserve 
an inequality since the denominator may be negative. However, because we are squaring the Euclidean 
norm, the common denominator is always a product of sums of squares. Taking the numerator is then 
equivalent to multiplying both sides by the denominator, a positive quantity, which will not change the 
direction of the inequality. Notice that the final implicant, ([7|), is simply Q, so proving Pq^x,k > for 
some K implies that x is GAS for the rational difference equation R. An algorithm for proving positivity 
will be shown in Section [3] Also note that whenever the function Q and equilibrium are clear from 
context, they will be omitted from the subscript of P. 

For a given Q and x we know that showing positivity of an associated polynomial implies GAS of x 
for R. We also know, given K, what that polynomial associated to Q and A' is. However, we still need 
to see how to conjecture a reasonable value for and then how to prove that the polynomial is indeed 
positive. We will see how to prove positivity in the next section. Now let's see how to conjecture a 
reasonable K value given R and x using a brute force method. Start with K = 1 and apply the following 
algorithm: 

1. Create the polynomial Pq^;y,k{'^) 

2. Apply a minimization technique to the polynomial Pq^;y,k{'^) (^-g-, simulated annealing, gradient 
descent, Metropolis-Hastings algorithm, etc.) many times to find approximate local minima of 
Pq,x,k- 

3. (a) If all approximate local minima are positive then conjecture that this K works, 
(b) If there is a negative minima then increment i^^ by 1 and go back to step 1. 

For ease of computation, and since this is only to conjecture a K, we apply the minimization technique 

in step [2] to a discrete set of points. We will restrict to a fine mesh with large upper bound. For example, 
fc+i 

the cartesian product X^_]^ {^5 ^e, . . . , Ne}, for some large value of N and small value of e. Then every 
point in the mesh is a vector of the form {iie, 22^, . • • , ifc+i^), where I < ij < N. 

Note that this is not the only possible algorithm for conjecturing a value for K. However, the main 
result in this paper is a positivity algorithm, so we will not consider other possible algorithms. One could, 
in theory, replace step 2 with the following, "Apply the polynomial positivity algorithm found in Section 
^\ Then step 3 would become, "// the algorithm in step 2 fails, increase K by 1 and go back to step 1, 
otherwise return K". Using this positivity algorithm, once a K value is found, it is also proved to be 
correct. However, using positivity in step 2 is sometimes not feasible since it often takes more computer 
memory than the conjecturing algorithm. 

3 An Algorithm to Prove Positivity of a Multivariate Poly- 
nomial 

So far, our algorithm to prove global asymptotic stability of a particular rational difference equation has 
reduced the problem to proving that an associated polynomial is positive. Now the question becomes, 
how does one prove positivity? In general one can show polynomial positivity using calculus, or using 
cylindrical algebraic decomposition [l]. However, both of these methods do not work particularly well 
when the polynomial has very high degree. This is typically the case for the polynomials produced in 
Section (2) so we must use a different algorithm. We propose a new algorithm which was inspired by the 
following definition and theorem found in 0. 

Definition 3.1. The polynomial P G M[xi, . . . , Xn] is 

• positive (resp. non-negative) from /i iff Vxi > /x, . . . , Xn > /i, P(xi, . . . , Xn) > (resp. P(xi, . . . , Xn) > 
0). 

• absolutely positive (resp. absolutely non-negative) from iff P is positive (resp. non-negative) from 
/i, and every partial derivative (of any order), P*, of P is non-negative from /i, i.e., Vxi > /i, X2 > 

/i, . . . , Xn > /i, P* (Xl , . . . , Xn) > 0. 

In addition, we will denote by cr^i,...,^^ (P) the polynomial obtained from P by translating in dimen- 
sion i by /ii in the negative direction. In other words, replace Xi by Xi + /i^ in P for all i. If /x^ = /x for 
all i then we simply write cr^(P). Also from [6 , a theorem that gives a necessary and sufficient condition 
for absolute positivity (and absolute non-negativity) is reproduced here. 
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Theorem 3.1 (Hong, Jakus 1998). Let P be a non-zero polynomial Then P is absolutely positive (resp. 
absolutely non-negative) from iff every coefficient in cr^{P) is positive, and the constant term is nonzero 
(resp. non-negative). In particular, if ii — then every coefficient in P is positive and the constant term 
is nonzero (resp. non-negative). 

Now, it is certainly too much to hope for the polynomials Pq^x,k be absolutely positive from zero. 
Of course, to satisfy Theorem it is only necessary that they be positive from zero (and possibly zero 
at a few points). Our algorithm will subdivide the positive orthant (the region in which all the variables 
are non- negative), denoted by R!j: where n is the number of variables in P, into regions in which P is 
positive on the boundary of the region, and essentially absolutely positive in some direction away from 
the boundary (i.e., there is a direction such that the directional derivative is positive). 

Since the polynomials we construct while pursuing global asymptotic stability are typically very 
complicated (high degree in all variables, some negative coefficients), we cannot easily show that a 
directional derivative is positive. Instead, for each region S C M+ we will create a polynomial, Ps{y) 
with the property that if Ps{y) > for all y then P{x) > for ah x e S. We will ffist describe the 

algorithm in two dimensions, and later generalize to the n-dimensional case. 

Let P := P{x,y) be a polynomial in two variables (n = 2). In order to show that P{x,y) > for 
(x^y) G M+, we first cut the positive quadrant into 4 regions as shown in Figure [l] where x is some 
positive number. In the case that P = Pq,x,k in Section [2] x will be the equilibrium point of the 
rational difference equation used to create P. 



y 



NW 



SW 



NE 



SE 



Figure 1: Cutting into 4 regions 

For each of these four regions we create a new polynomial from P by transforming the region into , 
and making the corresponding variable substitutions. See Figures [2]-[4]for the region transformations 
(transforming SE is analogous to transforming NW by permuting the x and y axes). Based on these 
region transformations we see that the associated polynomials are given by 

PNE{x,y) = cFx{P) = P{x + + x). 



Psw{x,y) = 

PNw{x,y) = 
PsE{x,y) 




The change of variables in P before applying a, inverting the variables or not, is self explanatory based on 
the transformation of the associated region. However, we must also multiply by and/or y^y (where 
dz = the degree of z in P ioi z = x,y) as needed before applying the a shift operator so that the resulting 
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X 



Figure 2: Transforming NE to 

Pn is still a polynomial. When talking generally about one of these polynomials, we will denote it by Pn, 
where the □ can refer to an arbitrary region. Note that if x = we will only consider the NE region, 
thus avoiding translating by | = ^ . 

Before we continue the algorithm by giving criteria to test positivity of the polynomials in ^ we 
must see why proving positivity of all Pn will be enough to prove positivity of P(x, y) for all (x,y) G M+. 

Proposition 3.2. Let P{x,y) be a polynomial, dx = deg^(P); dy = deg^(P); and x > 0. Consider the 
polynomials Pne, Psw, Pnw, Pse cls defined in If these four polynomials are all non-negative from 
then P(x, y) is non-negative from 0. 

Proof. For each of the four polynomials we will see that positivity for {x^y) G M+ implies positivity of 
P(x,y) in the corresponding region. 

If PNE{x,y) > for (x^y) G M+: Then by definition of PNE{x,y) we have 

Pne{x, y) = P(x + X, ^ + x) > for X > 0, ^ > 0. 

Let x' := X + X and y' \= y + x, then 

P(x^, y') > for x^ = X + X > X, and y' — y + x > x. 

This says precisely that P{x,y) > in the region NE. 
If Psw{x,y) > for {x,y) G M+: Again, by definition of Psw{x,y) 

for X > 0,2/ > 0. 

Following the previous case we first substitute x^ := x + 4 and y' y -\- ^ to get 

P f A , ^ J (x) "^"^ (y) > for x^ = X + i > i , and y = y + ^ > ^. 
\x^ 2/ / XX XX 

Since we are only interested in the region for which x^ and y^ are both strictly positive we may 
cancel the (x') " {y^ without reversing the inequality. We also make a second substitution letting 
x^ := ^ and y^^ := ^. Now we see that 

x' ^ y' 

Pix" ^y") > for < x^^ = ^ < X, and ^ < y" — — < x 

x' x' 

which is simply P(x,y) > in the region SW. 
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X 



Figure 3: Transforming SW to 



If PNw{x,y) > for {x,y) G M+: From the definition in ^ this means 

PNw{x,y) = p(^-^,y + x^ (^x+^^ >0 forx>0,2/>0. 

As in the SW case, we will make two substitutions. The first being := x + 4 and y' := 2/ + x. 
This gives us 

P f — , ^ ) (x) > for = X + 4 ^ 4 , and y^ — y -\- x > x. 

\x' J ^ ^ XX 

Again, we may cancel the (x^)^^ without reversing the inequality since x' is strictly positive in 
the region in question. Finally, we make our second substitution, x" := ^ (there is no second 
substitution for y') which yields 

P{x\ y') > for < x^^ = ^ < X, and y^ > x. 

x^ 

Therefore, P{x,y) > in the region NW. 

If PsE{x,y) >0 for (x,2/)eM+: This case is analogous to the NW case by interchanging the roles of 
X and y. 

In each of the four cases posit ivity of the polynomial corresponds to positivity of P(x, y) in the corre- 
sponding region. □ 
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X 



Figure 4: Transforming NW to 

To prove positivity of each of the Pn we will test two criteria, neither using anything more powerful 
than high school algebra. 

PosCoefFs: From Theorem |3.1| if all coefficients, including the constant term, of Pn are non-negative 
then Pu[x,y) > for {x,y) e M+. 

SubPoly: If the only negative coefficient in Pn (including the constant term) is on the xy term then we 
check whether the binary quadratic form, 

ax^ + bxy + cy^ (9) 

where a, 6, and c are coefficients of their respective terms in Pn, is positive definite (i.e., is positive 
for all {x,y) ^ 0) using its discriminant. The binary quadratic form discriminant of ([9| is defined 
to be d = 4ac — }? [12 . If a, d > 0, then (|9| is positive. Then, if this "sub-polynomial" of Pn is 
positive, Pn itself is positive (since the other coefficients are positive). Notice that this may not 
be the discriminant most are familiar with. A further discussion of why this is taken to be the 
discriminant can be found when the n-dimensional positivity algorithm is summarized later in this 
section. 

We also have an easy way to test whether P(x,y) < for some (x^y) G by checking the leading 
coefficient (the coefficient on the highest degree term) and constant term. 

LCoefF: The leading coefficient must be positive, otherwise the polynomial eventually tends to negative 
infinity in some direction. 

Const: Similarly, the constant term must be positive, otherwise the polynomial is negative in a neigh- 
borhood of the origin. 
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For each region □, if P\j passes one of PosCoefFs or SubPoly then, by Proposition 3.2 P(x,y) > 
in the region □. If Pn fails one of LCoefF or Const then we output false immediately because we know 
that there are points in region □ for which P(x,y) is negative. However, for some region □, if P\j has 
too many negative coefficients, its leading coefficient is positive, and its constant term is positive, then 
we must do more tests to establish positivity of Pn on . 

We would like to subdivide our original region (NE, NW, SE, or SW) into finitely many pieces and 
try again. However, there isn't an obvious way to do this since, except for SW, the regions are infinite, 
and we have used our "obvious" cutpoint, x. So instead, we will first map the infinite region into a finite 
rectangle with lower left corner at the origin (see figures |5] and [g]) and create a new polynomial Pn(x, y) 
from P{x,y) for each now finite region. We will then subdivide this finite region in order to prove that 
Pq(x, y) > on the region in which it is defined. These new polynomials will be defined in the following 
manner based on their corresponding region transformations. 





Figure 6: Transforming NW to finite rectangle (SE similar by interchanging axes) 



P^^(x,2/) = P(^,0x^^2/' 

PNw{x,y)^p(x,l^y^y 

P'sE{x,y) = P{ly) x^- 
P'sw{x,y) = P{x,y) 



restricted to < x, ^ < 4 

restricted toO<x<x,0<^< 4 
restricted toO<x< 4,0<^<x 
restricted toO<x<x,0<2/< 4 



(10) 



Along the lines of Proposition 3.2 we can guarantee positivity of P(x, y) given positivity of the related 
polynomials (10). 
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Proposition 3.3. Let P{x,y) be a polynomial, dx = deg^{P), dy = degy{P), and x > 0. Consider the 
polynomials P'ne^ P'sw ^ P'nw ^ P'se defined in (10). If any one of these polynomials, generally denoted 
P^, is positive on the region indicated in HOv, then P{x,y) is positive on the region □. For example, if 
P'nw{^^ y)>OonO<x<x and < y then P{x, y) > on the region NW , and similarly for the 

other polynomials/regions. 

Proof. We will only see the proof for P^y^(x, y), the rest will follow in nearly the same manner. Assume 
that P^v^(x, 2/) > for < X < X and < y < 4. Then by definition of Pj(^v^(x, y) we know that 

x^-^y"^"" >0 for < X < x,0 < 2/ < 4. 

yj X 

Since y is strictly positive in the region in which P^y^(x, y) is defined, we can cancel y^^ without reversing 



the inequality just as we did in the proof of Proposition 3.2 Now, let y' := ^ to see that 



P{x,y') >0 

for < X < X and y' = ^ > x, which is precisely the region NW. The other regions will follow by 
doing substitutions x :— ^ and := ^ as necessary. Note that no work needs to be done for SW since 
P'swi^^y) — P{x^y) cind is defined in the region SW. □ 

Next, we need to see how to prove that P^(x,y) > on the desired region. We will do this by 
subdividing the domain of into finitely many smaller rectangles. Then, for each smaller rectangle, 
S = {a < X < b,c < y < d}, we transform it to creating a corresponding polynomial, Ps{x,y), and 
test criteria PosCoeffs and SubPoly to see whether this polynomial is positive. See Figure [7| for the 
transformation of a general rectangle, to R^. Given this transformation, the polynomial, P^ , is given 
by 

P'six, y) = aii (p^ (^+a,-+c) x<y< 
b-a'd-c \ \x y J 



where Pg is one of Pne^ P'nw^ P'se^ P'sw^ ^^id d'z = degree of z in Pg for z = x,y. Before we see the 
canonical subdivision algorithm let us see why this Ps{x,y) will give the desired result. 

Proposition 3.4. Let P^{x,y) be a polynomial, d^ = deg^{P^), dy = deg^(PQ); < a < b, and 
< c < d. Consider the polynomial P's as defined in |TT] ). If Ps{x,y) is positive on then P^{x,y) 
is positive on the rectangle S = {a<x<b,c<y<d}. 



Proof. This proof follows the form of the proofs for Propositions |3.2| and |3.3| First, by definition of 
P'\S), the fact that Ps{x,y) > for x,^/ > means 

Pn ( / 1 + / 1 + c I {x+j^\ ^ f^+T^) ' ^0- 

V^+b^ y + ^c )\ V d-cj 

As in the previous proofs we may cancel ^x + -^zr^ ^ (^2/ + ^z^^ ^ without reversing the inequality. Let 
x' = x+ and 2/' = 2/ + ^ , then 

for x^ = x + -r^ > -r^ and y' — y-\- > -p—. Next, let x' = ^ and y'' = ^. Making this substitution 

o — a — b — a ^ ^ a — c — a — c ' x' ^ y' 

yields 

Pn {x" + a, y" + c) > 

for < x^^ = ^ < b — a and < y" — ^ < d — c. For the final substitution, let x" = x" + a and 
y'" = y" + c. Then 

7-,/ / /// ///\ \ r\ 

Pu[x ,y ) > 

for a < x'" = x" + a <b and c < y'" ^ y" + c < d. In other words, Pg (x, > for (x, y) e S. □ 
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d-c— 



b X 



b-t 




Figure 7: Transforming general rectangle, 5, to 



In principle any subdivision will work so long as we cover the entire finite rectangle. However, since 
the goal is to program the algorithm we need to specify a canonical subdivision. First we will simply 
divide into four equal regions. For each region we perform the above steps (transform the region and 
polynomial, apply criteria PosCoefFs and SubPoly). If we fail either criteria on a specific subregion, 
then we subdivide that subregion into four again and repeat. We continue to do this until we pass 
PosCoeffs or SubPoly, fail LCoefF or Const (and output false), or we reach some stopping condition 
and output FAIL. A stopping condition could be that we have subdivided N times, for some large N . 

Before we summarize the positivity algorithm in the general n-dimensional case we must see the 
general SubPoly criteria. It was previously stated only in the case of 2 variables. 

SubPoly-n: If the only negative coefficients in Pn (including the constant term) are on terms of the 
form XrXs then we check whether the quadratic form [10] , 



E 



i<j 



(12) 



where a 
all {xi 
defined as: A 



are coefficients of their respective terms in Pn, is positive definite (i.e., is positive for 
. . ,Xn) 7^ (0, . . . ,0)) using its corresponding matrix. The symmetric coefficient matrix is 



a. 



where 



■ a 



for all i ^ and ^ 



fT2). 



we can equivalently think of the quadratic form as xAx , where x = (xi, . 



Given this matrix, 
, Xn). Since A is a 



symmetric matrix, we know from the spectral theorem that it is diagonalizable by an orthonormal 
matrix, Q, so we have that A — QDQ^ where D is a diagonal matrix. We can then rewrite the 
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quadratic form as: 

From this we easily see that the quadratic form is positive definite iff all eigenvalues of A are positive 
(i.e., A is positive definite). Then, if this quadratic "sub-polynomial" of Pn is positive, Pn itself is 
positive (since the other coefficients are positive). 

We are now ready to summarize the algorithm in the n-dimensional case. Assume we have a polyno- 
mial P G M[xi, . . . ,Xn], and want to test whether P(xi, . . . ,Xn) > for (xi, . . . , Xn) G M+. 

1. First cut W\. into regions, similar to the NW ^ NE, SW, SE regions. For each variable we have 2 
possibilities for its domain 

< Xi < X or X < Xi < oo. 

A region is defined by making a choice for each variable, thus we have 2^ regions. The associated 
polynomial, Pn, for each region is then created by substituting 

( Xi -\- X ii X < Xi < oo 
x^hy^ ifO<x. <x 

in P, and then multiplying by (x^ + |)^^^ if < Xi < x. 

2. For each region we check our 4 criteria PosCoefFs, SubPoly-n, LCoefF, and Const. If all 2^ 
polynomials pass PosCoefFs or SubPoly-n then we are done, and P is positive on R!f:. If any of 
the polynomials fail LCoefF or Const then we are also done because we know that there are values 
in R!j: for the variables which make P negative. Otherwise, we continue on to step 3 for the regions 
which fail PosCoefFs and SubPoly-n. 

3. Assume we have a specific region (domains for each variable), R, which failed step [5] Then we 
create the polynomial P^ by substituting J- for Xi in P if Xi is restricted to x < x^ < cxd in i?, and 

then multiplying by x^"^* for those variables which were substituted. This new polynomial will be 
restricted to the region R\ which is defined from R in the following manner: if < x^ < x in i?, 
then Xi has the same restriction in R'; otherwise, x < x^ < cxd in i?, and then Xi is restricted to 
< Xi < 4 in R\ More formally, if D = {i : < Xi < X in i?}, and D = [n] ^ D then 



i?^ = ( X {0 < < x} 



x{o<..<i}) 



(a) Subdivide R' into 2^ equal regions, Sj, and for each region create the polynomial P's in the 



same manner as (11) 



(b) Test positivity of Ps- using criteria PosCoefF^ SubPoly-n, LCoefF, and Const If Ps. passes 
PosCoefFs or SubPoly-n then we are done in region Sj and can continue checking the rest of 
the subregions of R' . If P's^ fails LCoefF or Const then we stop altogether because we know 
that there are values for the variables in R!f: which make P negative. Otherwise, go back to 
step 3(a) with region R' now replaced by Sj. 

(c) If we have recursed more than N times (for some choice of N), stop and output FAIL. 

Before going on to the proof-of-concept for a specific difference equation and equilibrium let us see 
how to apply the polynomial positivity algorithm. We will see two examples, one in which SubPoly 
must be used, and one where subdivisions are necessary. 

Example 3.1. Let P(x,y) — x^ — xy + and x = 1. First we subdivide into the for regions NE^ 
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SW, NW, SE and get the following polynomials: 

PNE{x,y) ^ cri(P) = - + + X + 2/ + 1, 



Pswix,y) = cTi ( P ( i, i j x'^'^y'^y ] = - xy + y'^ + x + y + I, 

Pnw{x, y) = (71,1 ( P ( ^,2/ J a:"^"^ ) = x^y^ + 2x^2/ + 2x2/^ + x^ + 3x^ + ^/^ + x + 2/ + 1, 



Pse{x, y) - cri,i I^P l^x, - J J = x^r + 2x^2/ + 2xr + x^ + Sx?/ + 2/ + a: + 2/ + 1- 

The polynomials PArw(x, y) = Pse{x, y) have all positive coefficients, so they satisfy PosCoefFs. To see 
that PNE{x,y) (which equals Psw{x,y) in this example) is positive we must test criteria SubPoly. 

The only negative coefficient in PNE(x,y) is on the xy term, so we look at the sub-polynomial 
x^ — xy -\- y^ . The discriminant of this binary quadratic form is4-l-l — 1^ = 3, which is positive as 
needed. So we see that PNE{x,y) (and thus Psw{x,y)) is positive by SubPoly. 

In this example we don't have to do further subdivisions since Pnw, Pne, Pse, Psw are all positive. 
Therefore, we are done by Proposition |3.2| 

Example 3.2. Let P(x, y) = x^y — 5x^y + lOx^y + x + 2/ and x = 1. First we subdivide into the for 
regions NE, SW, NW, SE and get the following polynomials: 

Pne{x, y) — (Ji{P) — x'^y -\- x^ — x^y — x^ + x^y + 2x^ + 9xy + llx + 7^ + 8, 



1 1 



Pswix,y) = ai jx"^^/^ ) = x'^ + 4x^ + x^^ + 17x^ + 2x2/ + 21x + ^ + 8, 

PNw{x,y) - (71,1 [P [^,y] X 



x^y + x^ + Ax^y + 4x^ + IQx'^y + 17x^ + 19xy + 21x + 7?/ + ( 



PsE{x,y) = ai,i yP yx, - j y y j = X —x +xy + 2x + 2x2/ + llx + 2/ + 8. 

In this example we see that Psw{x,y) and PNw{x,y) pass criteria PosCoefFs since all coefficients are 
positive. For the other two regions we will need to subdivide because the negative coefficients are not on 
the term xy. 

Let us first examine SE. We need to create the polynomial P'sE^x^y) as in 
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PsE{x,y) = P ( -,y] X 



1 ^ 



— x^y + lOx^y + x^ — 5xy + y restricted toO<x< 1,0 < y < 1. 
Then we subdivide the region 0<x<l,0<2/<l into four equal rectangles: 

5i = |o<a;<^,0<y<^|, S2 = |o<x<^,^<2/<l|, 
l<x<l,0<y<^Y 54 = |^<a^<l,^<2/<l|, 
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and create four associated polynomials using (11): 

4 



P'^^ (x, y) + 3x^ + x'^y + + 4xy + 20x + % + 25, 

Ps^{x,y) =^x^y + 2^^^ + + + 3x^^ + 10x^+ 

25 

+ 10x2/ + 32x+ y2/ + 42, 

Ps^{x,y) =^x^y + Y^a:^ + 3x^2/ + 20x^ + 13x'^y+ 

+ 96x^ + 24x^ + 196x + 16^ + 144, 
Ps, (x, 2/) =§^^2/ + y + 10x^2/ + 34x' + 48x^y+ 

+ 166x^ + 98xy + 344x + 72^/ + 256. 

All four polynomials for the subdivision of SE are positive by PosCoefFs^ therefore P'se^^iV) ^ 0? cind 
by Proposition 3.3 we see that P{x,y) > on the region SE. 
Now let us look at Pne- Create Pne indicated by (10): 



P^E{x,y) = p(-^)x'-y'y 
\x yj 

— Sx'^y + + llx^y + 9x^ + 2x^y -\-x^ — xy — x-\-y-\-l 

restricted toO<x<l,0<^<l. 

Subdivide the region < x < 1,0 < y < 1 into the same 5*1, 5*2, S3, and S4 as above, and create the 
polynomials Ps., this time from P^e- 

P's-^ (x, y) —x'^y + 3x^ + 7x^y + 21^"^ + 19x^y + 58x^ + 
+ 33xy + 105x + 37y + 120, 



Ps,(x,y) 

Pssi^.y) 
PsA^,y) 



--^x^y + 4x^ + + 28x^ + 29x^y + 78x^+ 

105 

+ -Y^y + 144x + 6O2/ + 166, 

37 4 ^ 15 4^ 115 3 ^ 375 3^-,,. 2 _^ 
' 16^ 2/ + + y + + 142x 2/+ 

+ 463x^ + 320x2/ + 1040x + 272?/ + 880, 

15 4 ^ 83 4^ 375 3 ^.^n 3^ 463 2 ^ 
: — X 2/ + + -^x y + 130x + -^x y-\- 



+ 642x^ + 520x2/ + 1440x + 440?/ + 1216. 
As before, all four subdivision polynomials pass PosCoefFs, and so they are positive. Therefore, by 



Propositions 



3.2 



3.3 



and 



3.4 



we know that P{x,y) > for {x,y) G 



<> 



4 Proof of Concept 

We have now seen the full algorithm to prove GAS of equilibrium points of rational difference equations. 
However, there is no reason a priori that this algorithm is applicable. It could be the case that no such 
K (see Section [2] for a definition of K) exists, and this algorithm would be useless. 

We will now see that this technique does, in fact, work to prove global asymptotic stability of an 
equilibrium of a particular rational difference equation. The proof of the following theorem to establish 
global asymptotic stability will go through the procedure outlined in the previous sections. 

Theorem 4.1. For the rational difference equation 



4: + Xn 

1 + Xn-] 



(13) 



the equilibrium, x = 2, is GAS. 
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Proof. We will prove that K — satisfies Q from Theorem From the rational difference equation, 
the equilibrium x = 2, and = 5 we get the polynomial, P \— P{^2,2),5{{xi, X2)), as defined in (|6]): 

P = 25x?X2 + 340x1x1 + 1606x1x1 + 3060x?X2 + 2025xi + GOxiXs + 1158xiX2+ 
+8460xiX2 + 28936xiX2 + 45848xiX2 + 27090xi + TlxiXs + U18xtxl + 
+ 11229x1X2 + 53362xiX2 + 147345x1X2 + 207144xiX2 + 113103xi + 72xiX2 + 
+ 1420x1X2 + 9012x1X2 + 20174xiX2 + 24716x1X2 + 74718x1X2 + 163032x1X2+ 
+ 108952x1 + 47x^X2 + 1276x^X2 + \\V10x\x% + 25528x^X2 - 118780x^X2- 
-688300x^X2 - 1195361x^X2 - 790736x^X2 - 148969xt + 12x?X2 + 538x?X2+ 
+7854x?X2 + 45864x?X2 + 53604x?X2 - 515564x?X2 - 2066454x?X2- 
-2469564x?X2 - 207576x?X2 + 833882x? + x\x^2 + 86x^X3 + 2109x?X2+ 
+22070x?X2 + 102117x?X2 + 105526x?X2 - 695269x?X2 - 1867364x?X2 + 
+785343x?X2 + 6256056x?X2 + 4716817x? + 4xiX2° + 198xiX2 + 3530xiX2+ 
+29636x1X2 + 117218x1X2 + 136288xiX2 - 289440xiX2 + 253318xiX2+ 
+5674806x1X2 + 11634024xiX2 + 7054300xi + 4x2° + 148x2 + 2145x2 + 15348x2 + 
+53870x2 + 69340x2 + 30579x2 + 801874x2 + 3802411x2 + 6262908x2 + 3488704. 

The goal is to prove that this polynomial is positive when all variables are positive. Recall that we 
created this polynomial by taking the numerator of 

llAf- A'll' - ||Q'(Af) - A'll', 

where (5(A') is the map 





Xri 




r -1 


( 


)- 






Xn — 1 




Xn 



Now we run the polynomial positivity algorithm described in Section [3] to prove that this polynomial is 
positive. If the polynomial is positive when all variables are positive then the equilibrium, x = 2, is GAS 
for the original difference equation by Theorem |1.1| 

First we will prove that P > in the region NE. We make the polynomial Pne by substituting 
xi = xi +2 and X2 = X2 + 2 into P. See the Appendix in [5] for Pne and the rest of the polynomials 
as they will be omitted from this paper. Now we need to prove that Pne > in the region except 
when all variables are simultaneously zero. The only negative coefficient is on the term X1X2, so we can 
use the discriminant method. The binary quadratic form that we must show is positive definite is 

349366689X? - 6980904xiX2 + 318700575x2- 

Its discriminant is d = 445324725659927484 which is positive, so by SubPoly Pat^; > in R^. 

Now we will prove P > in the region NW. Create the polynomial Pnw by substituting xi = 1/xi, 
multiplying by Xi^^ = xf, and then translating xi by 1/2 to the left, and X2 by 2 to the left. All 
coefficients in Pnw are positive, and the constant term is zero. There is no proper subset of the variables 
for which setting them all equal zero yields the zero polynomial. Therefore, Pnw is zero only when all 
variables are zero, and so P > in NW. 

Next, we will prove P > in the region SE. First make the polynomial Pse by substituting 
X2 = l/x2, multiplying by X2^^ = X2°, and then translating xi by 2 to the left and X2 by 1/2 to the left. 
Now we need to prove that Pse > in the region R^ except when all variables are simultaneously zero. 
All coefficients are positive, and the constant term is zero. There is no proper subset of the variables 
for which setting them all equal zero yields the zero polynomial. Therefore, Pse is zero only when all 
variables are zero, and then P > in the region SE. 

Finally, we must prove P > in the region SW. Make the polynomial Psw by substituting xi = 1/xi 
and X2 — l/x2, multiplying by Xi^^ = xf and X2''^ = X2°, and then translating both variables by 1/2 
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to the left. Now we need to prove that Psw > in the region except when all variables are 
simultaneously zero. As in the region NE the term xiX2 has a negative coefficient (and that is the only 
such coefficient), so we will use the discriminant method again. The binary quadratic form that must be 
positive is 

349366689 2 _ 872613 318700575 2 

16384 ^04^ '''''' + 16384 

The discriminant is d = ^^^^^67108864^^^'^^ ' which is positive. Then, by SubPoly, Psw is positive in R+, 
so P > in the region SW. 

Since P > in all four regions, NE, NW, SE, and SW, the K value 5 is proven to work for the 
rational difference equation Xn+i = i^'l^'^ ^ ^ 

We can now see that the algorithm is indeed applicable. However, it wouldn't be possible without 
programming the algorithm. For large K values, even K > 3, the polynomials are near impossible to 
deal with by hand. For this reason there is a maple package, described below in Section [6] 



5 Results 

In this section we present the results that our algorithm can prove in full generality. So far we have con- 
sidered rational difference equations with specific numerical coefficients. In this section, we consider the 
case where the coefficients are additional variables, which are required to be positive. So the polynomial 
that we create is now a polynomial in the variables Xn,Xn-i, • • • ,Xn-k as well as all of the coefficient 
variables. We must point out that our algorithm will only apply when the equilibrium can be expressed 
as a rational function of the coefficient variables. For the proofs of the results found in this table see 5\ 
The equation numbers given in the following table match up with those in 2 , however the difference 
equations themselves may look different. The parameters presented here are to guarantee that the 
equilibria will be rational functions in the parameters. 



Eqn# 




Parameter Values 


Findings 


2 


^ 


M G M 


X = |M| is not LAS 


3 




M G M 


X = |M| is not LAS 


5 


/3Xn 


< /3 < 1 
1</^ 


X = is GAS 
X = is not LAS 


9 


'JXn-1 


< 7 < 1 
1 < 7 


X = is GAS 
X = is not LAS 


17 


1 M^-l 


M-1>0,M + 1>0 


x= |(M-1) is GAS 




4 l-\-Xn 


M-1<0,M + 1<0 


x = -|(M + l) is GAS 


23 


/3Xn 


< /3 < 1 


X = is GAS 






K/3 


X = /3 - 1 is GAS 


29 




< A < 1 


X = 1 — A is not LAS 




A+Xn 


1 < A 


X = is GAS 


30 




< A < 1 


X = 1 - A is GAS 






1 < A 


X = is GAS 


41 


a + /3Xn 


< /3 < 1 


X = ^ is GAS 


42 


^ ^ 4 Xn 


M -q<0,M + q<0,q>0 
M -q>0,M + q>0,q>0 


x = -^(M -q) is GAS 
X = |(M + ^) is GAS 


65 


1 M'^-q^+Axn 


M-q>0,M + q>0,q>-l 


x= ^{M -q) is GAS 




4 1+q+Xn 


M-q<0,M + q<0,q>-l 


X = -|(M + ^) is GAS 


109 




1 < A 


X = is GAS 


A-\-BXn+Xn-l 
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In addition to the results in the above table our algorithm can be used to prove GAS of many rational 
difference equations in which the coefficients have specific numerical values. See the Web Books on the 
web page that accompanies this paper 

http: //math.rutgers . edu/~eahogaii/GAS .html 

6 Maple Code 

In addition to the Web Books, there is a Maple package to accompany this paper which can also be found 
on the above web page. The three most useful procedures are ProveK, Prove, and WebBook. ProveK 
will use our algorithm to show that a given K value works to prove that the unique equilibrium of a 
particular rational difference equation is GAS. Prove utilizes ProveK to find the K value for a particular 
rational difference equation up to a given threshold. WebBook takes in a rational difference equation with 
variables for coefficients and proves GAS for a specified number of random choices for the coefficients. 
There is a Help function (type Help() to see a list of all procedures, and Help ((procedure name)) to get 
help on a specific procedure). In the help for each procedure, a sample is given for how to use it. 

7 Conclusion 

In Sections [2] and [3] we have seen both parts of our new GAS algorithm: first reducing the problem to 
proving that a polynomial is positive, and then proving polynomial positivity. Putting the two together 
we now have a completely algorithmic approach to proving GAS of a given rational difference equation. 
Inputs: 

R - rational function in A; + 1 variables 
X - equilibrium, solution to x = R{x, . . . , x) 
MaxK - a maximum K value to try 
Outputs: 

true if X is proven to be GAS for Xn+i = R{xn, Xn-i, • • • , Xn-k) 
false if X is not LAS for Xn+i = R{xn, Xn-i, • • • , Xn-k) 
FAIL if MaxK was not high enough. 
Procedure: 

1. Check local asymptotic stability using the linearized stability theorem. If not LAS then output 
false. If LAS then continue to Step [2] 

2. Conjecture a K value that satisfies Theorem 1 1 . 1 1 using the procedure outlined in Section [2] 

3. Apply the n-dimensional polynomial positivity algorithm outlined at the end of Section |3] If the 
conjectured K value was proven to work, output true. If the conjectured K value was proven not 
to work (Pk failed LCoeff or Const), or the algorithm reached a recursion limit, continue to Step 

HI 

4. U K < MaxK, increment by 1 and return to StepjS] U K > MaxK then output FAIL. 

This algorithm now gives us a completely automatic proof machine for global asymptotic stability. 
As was mentioned in the introduction, this problem has historically not been approached in any kind 
of systematic fashion. Many of the theorems found in [2] E] for proving GAS, were developed as gen- 
eralizations of techniques used to prove GAS of specific difference equations. This meant that given a 
particular difference equation, proving its equilibrium is GAS would amount to trying to apply various 
known theorems. Or, one may have to create a new theorem just to prove GAS of one particular rational 
difference equation. There may not have been a clear cut path leading to the proof. We believe that 
our new algorithm can serve as that path. Of course, given a difference equation that is known to be 
GAS, our algorithm may not always be able to prove it. However, we believe that it is much more widely 
applicable than any one previously known theorem guaranteeing global asymptotic stability. 
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